
% analog to NA_Dilation_multilevel_v4_2.m but for analysis of
% intra-cortical chrmine injected opto-stim experiments; basically defining
% the neurally active area as the mScarlet+ area in the widefield

% parameters for which data to use

vasomotionClass = 'maxDilation';        %'maxDilation','totalDilation','stimOnset', 'maxDilationTime' all options
validVesselSizes = [15,80];
usedistance = true;     %set this to true if you want to look at relationship with physical distance instead of distance weighted neural activity
explicitnegativecontrol = true;
modelgenotype = true;
correctForBaseline = true;            %set to true to use baseline value as a fixed effect in LME modeling
stimpower = 00;
plotMeanTrajectories = false;        % this will create figures with bootstrapped trajectories for each bin (not currently enabled)
numBootSamples = 1000;      %number of bootstrap samples to use; only relevant if you're plotting the mean trajectories, but if this number is high then function will be slow
plotcolor = 'm';
frameTime = .055;        %in seconds
prestimtime = 1.5;
% parameters for binning
numbins = 6;
minpoints = 5;
%% set which mice to use or not use
clear validmice;
%validmice = {'613'};
%validmice = {'981','932','929','904','933','621'};    %cre negative aEC cx knockout
%validmice = {'613','909','936','934','616','617'};         %cre positive aEC cx knockout
validmice = {'613','981','932','929','909','904','936','934','933','621','616','617'};      %all cx knockout mice
curpooled = aECcxOptoTable_histologyChrmine;          %will want to adjust this one...
curpooled.mouseID = categorical(cellstr(curpooled.mouseID));

baselineSubtract = true;        %this optional measurement will take the baseline value and subtract it from each trial. Mean baseline should already be zero, but since we're normalizing to the GLOBAL baseline, some trials have pre-dilated or pre-constricted vessels
if baselineSubtract
    intwindow = (27:134)*frameTime;
end
%% compile some of the relevant data:
validmice2 = curpooled.stimpower == stimpower;
if plotMeanTrajectories || baselineSubtract
    pooledtrajectories = totalIntensityTraj_normalized;
end

%% filter out datapoints not used in the current analysis:
VM = false(size(curpooled,1),1);
if iscell(validmice) 
    for n = 1:length(validmice)
        VM = VM | curpooled.mouseID == validmice{n};
    end
else
    VM = validmice;
end
if exist('validmice2','var')
    VM = VM & validmice2;
end
validmice = VM; clear VM
temp = curpooled.baselineWidth > validVesselSizes(1) & curpooled.baselineWidth < validVesselSizes(2);
validmice = validmice & temp;

if explicitnegativecontrol
    temp = curpooled.stimpower == 0;
    validmice = validmice & temp;
end
%% set evenly spaced bins to encompass all the data range in x:
if usedistance && ~explicitnegativecontrol
    temp = linspace(min(curpooled.dist2chrmine)-.01,max(curpooled.dist2chrmine)+.01,numbins+1);     %set bins based on the data range available
    temp = linspace(-.01,3000,numbins+1);           %set fixed bins
    numbins = numbins + 1;
    xbins = zeros(1,numbins+1);
    xbins([1,3:length(xbins)]) = temp;
    xbins(2) = 1;
elseif explicitnegativecontrol
    xbins = [-1,5000];
else
    xbins = linspace(min(x)-.01,max(x)+.01,numbins+1);
end
numbins = length(xbins)-1;
plotx = xbins(1:(length(xbins)-1))+(diff(xbins)/2);
y = nan(numbins,1);
ys = nan(numbins,2);
if modelgenotype
    genotypeeffect = nan(numbins,3);        %first column is the effect size, second column is the p value, third column is the confidence interval       
end

if baselineSubtract
    bl = mean(pooledtrajectories(:,3:26),2,'omitnan');
    switch vasomotionClass
        case 'maxDilation'
            curpooled.maxDilation = curpooled.maxDilation-bl;
        case 'totalDilation' 
            curpooled.totalDilation = curpooled.totalDilation- bl*range(intwindow);
    end
    pooledtrajectories = pooledtrajectories-repmat(bl,1,size(pooledtrajectories,2));
end
%% having recalculated the relevant values and filtered out the data you don't want analyzed, set a separate table for what you want
curpooled = curpooled(validmice,:);  
if plotMeanTrajectories
    pooledtrajectories = pooledtrajectories(validmice,:);
end
%switch resting diameter to the median value for that vessel across all the
%trials:
% get median resting diameter for each unique vessel:
[vid,~,ic] = unique(curpooled(:,[1,3]),'rows');           %MAKE SURE INDEXING IS CORRECT HERE!!!
for n = 1:size(vid,1)
    temp = curpooled.baselineWidth(ic == n);
    curpooled.baselineWidth(ic) = abs(median(temp,'omitnan'));
end
%mean correct resting diameter across ALL mice
curpooled.baselineWidth = curpooled.baselineWidth - mean(curpooled.baselineWidth,'omitnan');

%%  compute the average response in each bin:
% get bin-by-bin results
% can manually compute nested averages, variances or can use REML which is
% slightly more efficient but more opaque:
useREML = true;
alpha = .05;

for n = 2:length(xbins)
    ind = curpooled.dist2chrmine>xbins(n-1) & curpooled.dist2chrmine<xbins(n);           
    i2 = ~isnan(curpooled.vesselID);
    ind = ind&i2;
    if sum(ind) < minpoints
       continue;
    end
    i3 = ~isnan(curpooled.maxDilation);
    ind = ind&i3;
    curtable = curpooled(ind,:);        
    if size(curtable,1) < minpoints
        continue;
    end
    curtable.vesselID = categorical(curtable.vesselID);
    curtable.segID = categorical(curtable.segID);
    [~,~,ic] = unique(curtable(:,[1,11]),'rows');
    curtable.uniqueVessels = categorical(ic);       %this gives a unique number to each vessel such that different mice don't have overlapping IDs
    [~,~,ic] = unique(curtable(:,[1,3]),'rows');
    curtable.uniqueSegments = categorical(ic);
    
    if plotMeanTrajectories
        curtraj = pooledtrajectories(ind,:);
        
        IDcolumns = [1,13,14];      %this should be mouse - unique vessels - unique segments
        idx = nestedBootstrapNVCsample_3level_v2_0(curtable,numBootSamples,IDcolumns);
        numframes = size(curtraj,2);
        level1traj = zeros(numBootSamples,size(curtraj,2));
        for n2 = 1:numBootSamples
            level2 = ~isnan(idx(1,1,:,n2));
            level2 = sum(level2(:));
            level2traj = zeros(1,numframes);
            for n3 = 1:level2
                level3 = ~isnan(idx(1,:,n3,n2));
                level3 = sum(level3(:));
                level3traj = zeros(1,numframes);
                for n4 = 1:level3
                    temp3 = idx(:,n4,n3,n2);temp3 = temp3(~isnan(temp3));
                    level3traj = level3traj + mean(curtraj(temp3,:),'omitnan')/level3;
                end
                level2traj = level2traj + level3traj/level2;
            end
            level1traj(n2,:) = level2traj;
        end

        figure
        t  = (1:size(level1traj,2))*frameTime - prestimtime;
        cury = median(level1traj,1);
        curer = [cury-prctile(level1traj,2.5,1);prctile(level1traj,97.5,1)-cury];
        shadedErrorBar2(t,cury',curer',plotcolor);
        title([num2str(xbins(n-1)),' < range < ',num2str(xbins(n))])
        beautifyaxis(gcf);
        ylim([-6,70])
        xlim([-1.25,6.5])
        pause(.1)
    end

    % can build an empty model to calculate these parameters:
    if modelgenotype
        switch vasomotionClass
            case 'maxDilation'
                if correctForBaseline
                    lme = fitlme(curtable,'maxDilation ~ 1 + genotype + baselineWidth + (1|mouseID:movID)','FitMethod','REML');
                else
                    try
                        lme = fitlme(curtable,'maxDilation ~ 1 + genotype + (1|mouseID:movID)','FitMethod','REML');  
                    catch
                        lme = nan;
                        disp('not enough mice')
                        continue;
                    end
                end
            case 'totalDilation'
                if correctForBaseline
                    lme = fitlme(curtable,'totalDilation ~ 1 + genotype + baselineWidth + (1|mouseID:movID)','FitMethod','REML');
                else
                    lme = fitlme(curtable,'totalDilation ~ 1 + genotype + (1|mouseID:movID)','FitMethod','REML');
                end
            case 'stimOnset'
                lme = fitlme(curtable,'stimonset ~ 1 + genotype + (1|mouseID:movID)','FitMethod','REML');
            case 'maxDilationTime'
                lme = fitlme(curtable,'maxDilationTime ~ 1 + genotype + (1|mouseID:movID)','FitMethod','REML');
        end
    else
        switch vasomotionClass
            case 'maxDilation'
                if correctForBaseline
                    lme = fitlme(curtable,'maxDilation ~ 1 + baselineWidth + (1|mouseID:movID)','FitMethod','REML');          
                else
                    lme = fitlme(curtable,'maxDilation ~ 1 + (1|mouseID:movID)','FitMethod','REML');                 
                end
            case 'totalDilation'
                if correctForBaseline
                    lme = fitlme(curtable,'totalDilation ~ 1 + baselineWidth + (1|mouseID:movID)','FitMethod','REML');
                else
                    lme = fitlme(curtable,'totalDilation ~ 1 + (1|mouseID:movID)','FitMethod','REML');
                end
            case 'stimOnset'
                lme = fitlme(curtable,'stimonset ~ 1 + (1|mouseID:movID)','FitMethod','REML');
            case 'maxDilationTime'
                lme = fitlme(curtable,'maxDilationTime ~ 1 + (1|mouseID:movID)','FitMethod','REML');
        end
    end
    [~,~,STATS] = fixedEffects(lme);
    y(n-1) = STATS{1,2};
    feci = coefCI(lme,'Alpha',alpha);
    ys(n-1,1) = y(n-1) - feci(1);
    ys(n-1,2) = feci(1,2) - y(n-1);
    
    if modelgenotype
        if strcmp(STATS{2,1},'genotype')
            genotypeeffect(n-1,1) = STATS{2,2};
            genotypeeffect(n-1,2) = STATS{2,6};
            genotypeeffect(n-1,3) = abs(STATS{2,2}-STATS{2,7});
        elseif strcmp(STATS{3,1},'genotype')
            genotypeeffect(n-1,1) = STATS{3,2};
            genotypeeffect(n-1,2) = STATS{3,6};
            genotypeeffect(n-1,3) = abs(STATS{3,2}-STATS{3,7});
        end
    end
        
    
end

figure;
if explicitnegativecontrol
    errorbar(plotx,y,ys(1),[plotcolor,'s'],'markerfacecolor',plotcolor);
else
    shadedErrorBar2(plotx(1:6),y(1:6),ys((1:6),1),plotcolor);
end
    
if modelgenotype
    figure
    errorbar(plotx,-genotypeeffect(:,1),genotypeeffect(:,3),'ko-','linewidth',2,'markerfacecolor','auto')
    hold on
    plot([plotx(1),plotx(end)],[0,0],'k--')
    hold off
end
